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ABSTRACT 


Radar cross-section is a key element of low-observability. In order to 
reduce the cross-section of a particular platform, it may be necessary to determine 
the induced source distribution on the platform which produces the scattered 
electromagnetic radiation. Determining the distribution may be possible using a 
probe to measure fields on or near the outer surface of the object. However, the 
act of measuring may indeed influence the currents being measured. An alternate 
method is to back-propagate measurements made at distances beyond the realm of 
strong influence on the parameters of interest to construct visualizations of the 
local on-surface radiation contributions. This has been demonstrated for the case 
of cylindrical geometry. The theory is extended in this thesis to axisymmetric 


bodies for the special case of rotationally symmetric fields. 
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I. INTRODUCTION 


A. BACKGROUND 

1. The Importance of Stealth 

The concept of stealth has been around since long before man made his appearance 
on earth. Nature realized that the survival of its creatures depended on their blending in 
with the background. Evolution has provided the means for nature’s creatures to develop 
and refine stealth capabilities. It has been during the last fifty or so years that the military 
forces of the world have taken a keen interest in stealth [1]. The results were displayed 
by the success of the F117A stealth aircraft during the Persian Gulf War [2, 3, 4]. 

The success of these aircraft during the Gulf War can be attributed to many 
factors. By far the most important factor was possession and use of stealth technology. 
The U.S. military had it and the enemy did not. More importantly, the enemy did not 
possess the technology to defeat a stealthy airborne platform. That was in 1991, more 
than six years ago. Since then stealth technology has been ever so slowly let out of its 
“black box.” It is no longer “our own little secret weapon” that only the U.S. military 
possesses, and there is no guarantee that during the next conflict our adversaries will not 
possess anti-stealth technology. In fact ultra-wideband radar (UWB) was being 
developed for detection of stealth platforms before the Gulf War began [5]. In 1994 the 
United States Air Force admitted that some radars, including some mobile units, could 
detect the B-2 bomber [6]. And as one would expect, the Russians are also developing 
their own anti-stealth technology [6]. The bottom line is that if stealth platforms are to 


remain stealthy, they need to continue to decrease their observability. 


The idea of stealth is comprised of several aspects of observability including 
infrared radiation, optical, acoustic, and radar echo. The objective is to minimize the 
observables associated with your platform so as to reduce the chances of being detected 
by the enemy. While current stealth platforms have incorporated reductions in each of 
these areas, the most emphasis has been placed on reduction of radar signature. This is a 
consequence of radar being the primary means of detection for the U.S. military, its allies, 
and its adversaries. That being the case, this thesis attempts to take one small step in the 
direction of radar cross section reduction by determining the induced source distribution 
on a body which produces scattered electromagnetic radiation. 

2. Goals of the Research 

Radar cross-section (RCS) reduction is a key element of low-observability. 

Current techniques used in determining the RCS of a platform rely on analysis of 
measurements made in the far-field. Generally these measurements provide a gross picture 
of the platform’s overall RCS as a function of viewing angle. This information enables 
engineers to then modify the design in an attempt to reduce its RCS. This becomes an 
iterative process with design modifications leading to measurements which lead to more 
design modifications. 

At the other end of the spectrum, measurements made using a probe on or near the 
outer surface of the body may influence the quantities being measured. In other words the 
act of measuring would induce errors into the quantities being sought. 

The focus of this study is to more accurately determine the source of scattering on 
a body. This will be done by evaluating the viability of back-propagating measurements 


to an axisymmetric body to image the source of scattering on the body. This analysis will 


enable more precise location of scattering sources than does gross analysis of the far-field. 
Once the local scatterers have been identified it is then theoretically possible to remove or 
reduce them in the overall effort of radar cross-section reduction. 

The objective of this study is to investigate the viability of back-propagating 
electromagnetic field measurements to an axisymmetric body in order to determine the 
source distribution for radiated power from the body. If this is indeed possible, then the 
follow-on objectives include determining the range and spatial resolution at which 


measurements must be made in order to provide meaningful results. 


B. SOURCE IMAGING 

1. Acoustics Work with Cylinders 

A brief history of this subject starts with researchers attempting to determine 
sources of acoustic noise on a body. It was shown that acoustic intensity could be used to 
localize sources of sound on structures which radiate to the far field. A new quantity 
called the supersonic acoustic intensity vector was defined and its application to 
measurements on plate and cylinder-like structures was demonstrated [7]. 


“Supersonic intensity is composed only of wave components which 
radiate to the far field (supersonic), with the non-radiating (subsonic) 
components eliminated. The normal component of this supersonic intensity 
vector, measured in the extreme near field or on the surface of the 
structure, provides an accurate tool for locating regions (“hot spots”) on 
the structure which radiate to the far field. Furthermore, the supersonic 
intensity provides an accurate quantification of these source regions, 
providing a ranking of the strength of the identified source regions as a 
function of frequency. This identification and ranking provides a powerful 
new tool in the understanding and control of radiated noise.” [7] 


Coupled with the finding that acoustic source regions could be determined through 
back-propagation was that the results were obtained in two dimensions. Using a cylinder 
as the scattering body, scattering data measurements were made on an arc of a circle with 
diameter on the axis of the cylinder [8]. In the spherical (7,0,¢) coordinate system this 
translates to measurements over 9 at aconstant ¢. The restricted case of axisymmetric 
fields and currents is the next level of difficulty before working with full three dimensional 
shapes. 

je Extension to Electromagnetics for Cylinder Geometry Case 

The discoveries employing supersonic modes to enhance acoustic imaging on 
cylindrical shells has been translated into the realm of electromagnetics for the case of 
cylindrical geometry. It was shown that superluminal modes which satisfy |k,| < k and 
have faster-than-light axial propagation velocity provide time-average power to the far 
field. An optimal propagation deconvolution filter was implemented to remove the effects 
of the subluminal modes during back propagation. [9] 

SF Extension to General Non-separable Geometries Using Finite Element 

Methods 

The goal of this research is to extend the concepts developed in previously 
conducted research to bodies with non-separable geometries. By definition these bodies 
have shapes that do not lend themselves to explicit mgorous solutions using separable 
modes, such as cylindrical or spherical wave functions. A finite element method (FEM) ts 
used to relate the conditions at the boundary (measuring region) to conditions on the 


surface of the body, in this case relating the scattered field to the source current that 


generated it. The number of elements used is chosen so that the contour of the body is 
closely matched. Figure 1 shows a meridian plane of a sphere surrounded by a finite 
element mesh which connects the surface of the sphere to the arc of a circle with diameter 
on the axis of the sphere, as previously mentioned. In Figure 2, the body represented is a 
cone. The more contours on the body, the greater the number of elements required to 
closely match those contours. Ultimately the accuracy of the solution is determined by 


the number of elements used. 


2-D Axisymmetric Shape with mesh 


Mesh Parameters: 
Nodes: 
lements: 544 
nknowns: 264 


meters 





0 2 4 6 
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Figure 1. Sphere (Meridian Plane) Surrounded by Finite Element Mesh 


2-D Axisymmetric Shape with mesh 


Mesh Parameters: 
Nodes: 
Elements: 5/8 
Unknowns: 282 


meters 





meters 


Figure 2. Cone (Meridian Plane) Surrounded by Finite Element Mesh 


In the chapters that follow, the source current on the surface of an axisymmetric 
body is determined through back-propagation of the scattered field it generates. In 
Chapter II, the scattered field (measured field) is determined through numerical integration 
of a surface current induced - an axisymmetric body by a locally placed elemental dipole. 
Chapter III details the back-propagation of the scattered field. Chapter IV provides an 
analysis of both the numerical integration results and the back-propagation results. 
Conclusions are presented in Chapter V. An appendix which includes all pertinent 


computer code is also included to assist in the understanding and continuation of this 


research. 


IH. GENERATION OF SCATTERED FIELD (MEASURED FIELD) 


A. BACKGROUND 

Before the scattered field can be back-propagated it must be measured. To 
simulate the process and control the sources of error, the measured field is computed to a 
specified level of accuracy. While reducing potential error associated with field 
measurements, generating fields leads to an additional hurdle, namely how to generate a 
scattered field from an arbitrary axisymmetric body. This problem is solved by inducing a 
surface current on the body and integrating it to find the resulting scattered field. 
Potential error is introduced as a result of integration accuracy. A sphere is first chosen as 
the arbitrary axisymmetric body in order to test the accuracy of the general source 
integration algorithm. 

A particular scattered field is generated by placing a radial directed elemental 
dipole in the vicinity of the metal sphere. The dipole induces a surface current a on the 
sphere which in turn generates scattered fields E. and H. _ The measured field H at 
some specified distance from the sphere is composed of the field scattered (H, ) from the 
sphere, and the field incident (H,) from the elemental dipole. The ‘exact’ induced 
currents and scattered fields can be computed for this test case and used to validate the 


accuracy of the numerical integration and finite element algorithms. Calculation of the H 
field is considered in Section C. It is generated using the program DIPSPHR2.M found in 
the Appendix. Field integrations for axisymmetric currents are detailed in Section D and 


integration results are provided in Section E of this chapter. 


B. SPECIALIZATION TO AXISYMMETRIC FIELD CASES 

Determining the source distribution on bodies with non-separable geometries is the 
ultimate objective of research in this area. This thesis reduces the problem from three 
dimensions to two dimensions by assuming axisymmetric fields generated by axisymmetric 


currents on a body of revolution. There are two special cases to be considered: the TE, 


case in which the E field and the surface current are transverse to d , and its dual, the 
TM, case in which the H field is transverse to d and the surface current is d -directed. 


These cases are depicted in Figure 3. 


E, 
= H,>) =5 yee 


Zz 


Restricted Case: 7E, Dual: 7M, 
Figure 3. Axisymmetric Field Cases 


This thesis will investigate the 7E, case. The fields for this special case are, 
H(p,z)= H,(p,2)¢ (1a) 


E(p,z) = E,(p,z)p+E,(p,z)2 (1b) 


The fields for the dual special case, 7M,, are, 
E(p,2) = E,(p,2)¢ (2a) 
H(p,2) = H,(p,2)) + H,(p,2)2 (2b) 
Maxwell’s curl equations are: (with @ = 27 ) 
VxH = joceE+J (3a) 
VxE=-jopH (3b) 


Using Equation (3a) for the assumed field components in Equation (1) gives 


_~la, ay 4 
a a a 2 a) 
OE seh, ibis Ab 
pn ae) ee (4b) 


These are the generating equations for E interms of H , and known source current on 


Using Equation (3b) for the assumed fields in Equation (2) gives 


: ¢ 

LS oe (Sa) 
i= ae (5b) 
On. = 

H upp” 


These are the generating equations for H interms of E a 


Cc OFFSET DIPOLE TO GENERATE SURFACE CURRENT ON SPHERE 
A radial directed dipole located in the vicinity of a metallic sphere is used to 
generate surface currents on the sphere. These surface currents are then integrated to find 


the scattered field at a specified distance. The dipole case is also used to test the 


integration accuracy by providing a nearly exact scattered field. Figure 4 shows an 
elemental dipole located in the vicinity of a metal sphere in the spherical coordinate 


system. 








z-Directed Elemental 
Dipole at (0,0,z,) 


Metallic Sphere of 
radius @ __ we 





Figure 4. Dipole in Vicinity of Metal Sphere 


The procedure for determining the surface current consists of several steps. The 
E. and H, fields due to the dipole are represented in (7,9,9) coordinates of the centered 
sphere. The scattered fields E. and H. due to the induced currents on the sphere are 


then expanded using spherical harmonic expansions with unknown coefficients. Enforcing 


?x(E,+E)=0 => Ewan = 0 on the sphere surface allows solution for the expansion 
coefficients. Finally ve =a (H, +H.) is used to find induced surface currents which 


produce E. and H, fields scattered from the sphere. 


An elemental dipole placed at the origin produces the following fields [10]: 

















7] Idl 1 J — jJPry 
E, = a 00884 7 Br: I JB (6a) 
Noldl . IB J TT 
= —— 4 6b 
Eg, 4 810 0.) _ “ aa e (6b) 
Idl 1 
3 47 Pe if 


The position of the dipole must then be translated to the position z = z,. This 


translation is accomplished using the following transformation formulas: 


r, =r? +z) — 2rz, cosO (7a) 





A= 
cos@, = id (7b) 
rT, 
in@ 
sin@, = = (7c) 
d 
E, = £,, cos(@ — 6,) — E,, sin(@ — @; ) (8a) 
E, = E, sin(@- 0,)+ E,, cos(@ — @,) (8b) 
ot (8c) 


Pl 


The transformation formulas in Equations (7) and (8) are then used with 


Equation (6) to express the fields due to the dipole on the surface of the metal sphere, 
E,;(a,@), E,,(a,@) and H,;(a,@). This field will induce a surface current J = J,(0)0 


on the surface of the sphere which will, in turn, generate scattered fields E,., E,,, and 


r 


H,,. The tangential E field at the sphere surface will exactly cancel the tangential E 


field due to the dipole thus producing a total tangential E field of zero 
FE, ,(a,0) = —E,; (4,9) (9a) 
The surface current will also satisfy J =" Hoa On the surface: 


J =? x (Hy, +H,,)0|_ =- (Hy, +H,,)6|_ xp) =- 


=> Jq(8) = -(H,,(a,0) + H,,(a,8)) (9b) 
This J,(@) is what needs to be integrated to find the radiated H, (7,0) at locations off 
the surface of the metal sphere. The result is H,, which will be used to compare results of 


the integration covered in the following section. 


Note that although J, is produced by the total H, on the metal surface, it 
generates only the scattered field. The scattered field has an E,, which exactly cancels 
out the E,, of the incident dipole field on the sphere surface. 


To determine the “exact” solution, the scattered fields can be expressed as 


weighted sums of spherical harmonics [11], 


«0 A 
iE U.o) — ith 24, — An Pr) ps '(cos@) (10a) 


2 


wo (2) 
H,,(r,0) = da, ——" (cos6) (10b) 


where H® (Br) = if (Br) - ne (Gr) is the “Riccati” spherical Hankel function of the 
second kind representing outbound waves. P'(cos@) is the “associated Legendre 
function” and is the m=1 case of P”(cos@). To numerically solve the problem, 


truncate the series in Equation (10) and solve for WN values of complex a,. Using 


He” 
Equation (9a) and defining D, = — gives 
N 
J 2,9,D,P, (cos 6) = -E,,(a, 8) (11) 
n=] 


The a, are found by use of orthogonality of the associated Legendre functions, 


P™(cos@). For the special case m= 1, 


n 2n(n + 1) ] = 
| P'(cos6)P!(cos0)sin@dO=4 2n+1 (12) 
0 0 ln 
We can sift out coefficients: 
FDI) = — | E,,(a,0)P! (cos6) sin dO (13a) 
2n+1 4 
== = esa ae 
I 2n+1 
__ fn (antl) (13b) 


acs 
" JmD, 2n(n+1) 
The /, integration can be performed numerically to any accuracy since /,,(a,@) and 


P'(cos@) can be obtained using as many @-points as desired. Once the a, are found the 


13 


J,(@) can be computed at any @ value by use of Equation (9b) with H,, given by 


Equation (1 0b). 


D. FIELD INTEGRATION FOR AXISYMMETRIC CURRENTS 

As noted in the introduction, the first step is to measure the field that will be back- 
propagated to the surface of the axisymmetric body. For our purposes we will assume an 
axisymmetric surface current distribution on the body. The “measured” fields will actually 
be determined through the integration of the surface current. In this context the scattered 
field is generated through integration rather than through illumination and measurement. 

The surface of an arbitrary axisymmetric body of revolution can be defined by use 
of the generating contour p’(z’) using the (9,¢,z) cylindrical coordinate system as 


shown in Figure 5. Source points on the body are defined by individual (p,z) pairs. The 


field point in space is located at (p,9,z). The axisymmetric surface current distribution is 
J(r') = J (p’,2’t +J (9'.2')6 . The unit vectors ¢ and gd are tangent to the surface at 


each point and both J, and J,are constant with changing ¢. 
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surface 


Figure 5. Axisymmetric Body of Revolution 


The fields will be numerically evaluated via the vector potential formulation given 








by [11]: 

aie l oe 

H(r)=—VxA(r)_. (14a) 

by 
a efit 1 wy A) 14b 
a /' aiieaiialaaeaaar TE (14b) 

where 

== Me tps me tes 

A(r) = az ee) = as (15) 
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with R=r-r’ and R= Ir _ r'| . For field points and source points with coordinates 
(p,¢,z) and (p’,¢’,z’) the law of cosines gives 

R= p" +p” -2pp' cof d - 9’) +(z-2') (16) 
The curl operator in Equation (14a) 1s taken inside the integral in Equation (15) and note 


; — JPR : ee Pe centers 
that it only operates on the © Up term (unprimed position in R = r = r’| ) 





HG) =— [|v [= ; “| Fas (17a) 


WC are 


where 





(Set) = ae - 


Thus the general solution is: 


FE) a — [7 (r') « RY LPR) eras: (18) 


surface 


The next step is to derive the explicit integrations to be performed. To do so we 
select 7 in the x-z plane, with ¢ = 0as shown in Figure 6. Since oe and the resulting 
fields so derived are axisymmetric (not functions of ¢ ) we lose no generality by doing this. 
In fact, if we find H,,H, and H, at this field point we note that for other locations where 
g # 0 the substitutions H, > H,;H, > H, and H, > H, apply. Also note that 

r= pxr+22 (19a) 


r’= p'cos@’x + p'sing'y + 2'z (19b) 
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R=r-r'=(p-p'cos¢')£- p' sing’) +(2-2')2 (19c) 


R=|Rl= Vp’ +p” —2pp' cos¢’ +(z-2z')’ (19d) 


N > 
\ 


Field Point 
(x,0, z) o 


ae (p', 9',2’) 





Figure 6. Geometry for Explicit Integrations 
The field integrations will be performed numerically by breaking the axisymmetric 
surface into a large number of flat circular rings as shown in Figure 7a. These rings are 


stacked to approximate the surface of revolution. Each ring 1s “flat” in the direction of 


t and curved in the direction of d as depicted in Figure 7b. 


i 


Stacked Circular 
Rings 


= A 


Figure 7a: Stacked Circular Rings form Axisymmetric Body 
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“Flat” in ¢ direction 
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Figure 7b: Single Circular Ring 


Figure 7. Decomposition of Axisymmetric Body 
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Consider now evaluating H(r) in Equation (18) due only to J(r’) = J,(s)t over 
a single ring as in Figure 7b. Over the ring we will approximate /,(s) as piecewise linear, 


specified using its values at the top and bottom of the ring, J,(s,) and J,(s,) via 


(s—s,) (s-s,) 


va ae na eds —s,} (20) 








This equation 1s a straight line interpolating function between known values at s, ands,. 


A single ring is shown in Figure 8a where s 1s the path length over the ringed surface from 


the + z axis. The piecewise linear variation of /,(s ) between s, and s, is shown in 
Figure 8b. 
To evaluate Equation (18) for J(r’) =J, (s)t over a single ring, we note that 
t = cosacos¢’X + cosasing’ y —sinaz (21) 
where cosa =/-p and —sina=/-z. After simplifying, 
ix R=9t,R, -1,R,)+3(t,R, —1,R,) + (t,R, —,R,] (22) 


and substituting Equation (22) into Equation (18) with restriction to one ring gives, 


HL = a Je — z')cosa@ — p’sin a|F,(s)p'ds (23a) 
H, = ra we (s)[p" sina@—(z—-2’') cosa|F, (s')—psinaF, (s)}!p'ds (23b) 


a: 
H, = —|J,(s)pcosaF,(s)p'ds (23c) 
4 : 


s 3O 






Single Ring 


Figure 8a: Single Ring 


J,(s) 


J,(s2) 


Ji(s;) 


wf 
— 
172) 
No 


Figure 8b: Piecewise Linear Variation _—_J,(s) 


Figure 8. Surface Current Variation Over One Ring 








where 
os 1+ jBR) _, 
F,(s) = f sing z Je Mag’ (24a) 
{ 1+ jBR\ _, 
P(s)= fos R Je PR dg! (24b) 
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F(s)= | (2 1eJPR) imag’ (24¢) 


As we vary 9’ while holding s constant (e.g. p’,z’ constant), R(¢’)is an even 


sw 





function: R(—¢’) = R(g’). Therefore (Ht Je “J +s also an even function of ¢’. 


Since sing’ is an odd function of 9’ the integrand in Equation (24a) is also odd giving 


F,(s) = 0. Weare left with only H,, in (10b) since H, = H, = 0. 


Extending this to any position r, we have for a SG 


H=H,(r)¢ (25a) 
where 
H, = = [ 4,(9){[¢' sina - (2 - 2") cosa] F,(s)- psinaF,(s)}p'ds 25») 
and 
F,(s)=2 cos¢' (tt a tI, sR ag! (26a) 
eee (r= v JPR “iB ag! (26b) 


E. CENTERED SPHERE INTEGRATION TEST CASES 

A centered sphere was chosen as the test case to determine the accuracy of the 
integration algorithm. The integration algorithm was developed using the trapezoidal rule. 
Nine cases were tested with sphere sizes ranging from a radius of one wavelength to ten 
wavelengths. For each sphere size the H field was generated at three different locations: 


0.2, 1 and 3 wavelengths from the outer surface of the sphere. 


Zl 


The offset distance of the dipole was kept constant for each of the three cases associated 
with the three different sphere sizes. The number of modes used was determined by the 
following: 

N=integer(2 ka + 2), (27) 
where & is the wave number (27/A) and a is the sphere radius. In each of the nine 


cases the frequency of operation was 300MHz. The RMS error was determined by 


comparing the H field determined through integration to the H field exact solution found 
using the procedure developed in Section II.C. Each of the cases was tested with 
increasingly fine segmentation in @ and ¢ on the sphere surface until the RMS error was 
in the neighborhood of one percent. In order to achieve the desired integration accuracy, 
the number of integration points on the body was increased. This relates to the spatial 
resolution required when measuring the H field in order to achieve acceptable results 
when back propagating that field. As the number of integration points on the body 
increases, the measured H field becomes more accurate. The more accurate the 
measured field, the better the surface current on the sphere can be resolved. The test 
cases summarized in Table 1 are chosen based solely on the RMS error associated with 
each being in the neighborhood of one percent. 


The results shown in Table 1 indicate the RMS error for source integrated 17, 


observed at the “field distance” radius. The RMS error is a function of the number of 
surface points at which the surface current was integrated. Figures 9 through 12 show 


how the integration results converge for several of the cases shown in Table 1. For the 
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CENTERED SPHERE INTEGRATION RESULTS 


Case | Sphere | Field Dipole | Modes 
Radius | Distance] Segments | Segments | Offset RMS Error 
(7%) 


(A) (A) 
1.1490 
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=p 
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1.1620 


1.1510 


1.0310 


Ws 


0.1220 


0.1250 


0.8268 


pt 
© 
N 





0.0456 


0.0459 
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fused _—" 
lod — 





Table 1. Centered Sphere Integration Results 


cases in Table 1, the number of ¢ and 6 segments is based on the radius of the sphere. 
However, the number of ¢ segments remains constant as the radius of each ring (see 
Figure 7) decreases near both poles of the sphere. This finer ¢ segmentation gives 


increasingly narrower patches of surface area over which the integration is performed. 
Consequently, the RMS error does not appear to decrease significantly as the number of 


g segments is increased. This is a result of the increased ¢ segmentation that is “built in” 
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as the radius of the rings decrease. An increase in RMS error is observed for the 


decreased field distance from the sphere in Cases 2a and 3a. Since the field distance is 


only 0.2 from the sphere surface in these cases, the large relative contribution and rapid 


variation of inverse distance terms in the integrand of Equation (18) for surface points 


closest to the field point require increasingly fine segmentation for accurate numerical 


integration. 
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Figure 9. Integration Convergence: Casela 
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Integration Convergence: Case1b 
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Figure 11. Integration Convergence: Case 2a 
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Figure 12. Integration Convergence: Case 2b 


The magnitude and phase of the H field at the measured distance are presented in 


Figures 13 through 30. These figures help to tell the complete story of the computed H 
field. In each figure the “exact” quantity computed via series solution is depicted by a line 
and the “measured” quantity is depicted by dots. Note that in both the magnitude and the 


phase plots the measured quantity closely matches the exact quantity. These results 


indicate that the integration produces an H field with minimal error. This knowledge will 
help isolate the source of error that may present itself due to the back-propagation 


process, which will be described 1n Chapter HI. 
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Figure 13. Magnitude Comparison: Case la 
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Figure 14. Phase Comparison: Case la 
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Figure 15. Magnitude Comparison: Case 1b 
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Figure 16. Phase Comparison: Case 1b 
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Figure 17. Magnitude Comparison: Case Ic 
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Figure 18. Phase Comparison: Case lc 


Hexact (line) & Hinteg (dots): Case2a; RMS Error= 1.031% 
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Figure 19. Magnitude Comparison: Case 2a 
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Figure 20. Phase Comparison: Case 2a 
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Figure 21. Magnitude Comparison: Case 2b 
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Figure 22. Phase Comparison: Case 2b 
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Figure 23. Magnitude Comparison: Case 2c 
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Figure 24. Phase Comparison: Case 2c 


Hexact (line) & Hinteg (dots): Case3a; RMS Error= 0.8268% 
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Figure 25. Magnitude Comparison: Case 3a 


Hexact (line) & Hinteg (dots): Case3a 


i 


2) 





0 50 100 T50 200 
Theta (degrees) 


Figure 26. Phase Comparison: Case 3a 


33 


Hexact (line) & Hinteg (dots): Case3b; RMS Error= 0.04562% 
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Figure 27. Magnitude Comparison: Case 3b 
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Figure 28. Phase Comparison: Case 3b 
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Figure 29. Magnitude Comparison: Case 3c 
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Figure 30. Phase Comparison: Case 3c 
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HI. BACK PROPAGATION OF SCATTERED FIELD (MEASURED FIELD) 


A. BACKGROUND 

A finite element method is used to back-propagate the measured field to the 
surface of the sphere. The measured fields that are back-propagated are those determined 
in Chapter II and summarized in Table 1. Once the measured field is back-propagated, the 


surface current is determined using Equation (9b). The surface current due to back- 
propagation OF ax: ) is compared to the original surface current To. ) due to the 


elemental dipole which is computed using a spherical harmonic series, per Section II.C. 


The following section describes the finite element formulation. 


B. FINITE ELEMENT FORMULATION FOR AXISYMMETRIC FIELDS 


For the special case being considered, with no @¢ variation in materials or sources, 


and J, =0, the fields are derived in Section II.B., results of which are repeated here for 
convienence, 
H(p,2) = H,(p,2)¢ (28a) 
E(p,2) = E,(p,2)p + E,(p,2)2 (28b) 
Maxwell’s curl equations are: (with @ = 277 ) 
VxH= jocE+J (29a) 


VxE=-jopH (29b) 


S7 


Using Equation (29a) for assumed fields in Equation (28) gives 


oe Bo ap (30a) 
, 110 1 
J@ a pe ee (30b) 


These are the generating equations for E in terms of H , and known source 


current J. In this case Equation (29b) simplifies to: 
ee. . ere 
JO(V x E)-¢ = jo : es =@ *yH, (31) 


Substituting the E fields from Equation (30) into Equation (31) gives the partial 
differential equation (PDE) satisfied by H,: 

4 (oH vip + 2 1 om, = SC) CI,) G2) 
where €= €(p,z) is allowed to be to be inhomogeneous. 

The PDE in Equation (32) will be ee solved using the finite element 
method (FEM), where H, is approximated with electrically small triangular element 
regions in the (p,z) plane. [12] 


A variational Euler-Lagrange algorithm is used to implement the finite element 
method. This algorithm replaces the solution of the PDE in Equation (32) with finding the 


stationary point in complex function space of a quadratic (in H, ) functional. We first 


expand the field using pyramidal basis functions, u,(p,z), defined for each node in the 


mesh and having support of the surrounding tnangular elements, 
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H,(p,2) = Lihrty(P,2) (33) 


and then substitute Equation (33) into the functional. The stationary solution is found by 


differentiating with respect to the unknown values of H,, located at the m-th nodes, to 


yleld the linear system 
= ] 
ar, Sf —V(eu,)-V(pu,) - Ppu,u,dpde = 0 (34) 
n=] overlapping elements p 


Defining the double integrals indexed on m,n as T(m,n) , the linear system in Equation 


(34) can be rewritten as, 


> hT(mn)=- > h Timn) (35) 


n inside boundary n onboundary 


—na(n) hb(n) 


0] h, 
“| A h, || Bye G6) 
0 +d Ay 


St m 
Ny * Nyy Square Array Nu xNs 
(very sparse) (sparse) 


The m-th row of A and B is formed by Equation (35) with the n-th term 
corresponding to a column of A if h, isunknown. The appropriate column of B is filled 
if h, is a boundary node. There is no contribution to A or B for nodes on the z-axis 


(0 =0) since H,(p=0)=0 and those nodal values of h, are known. [13] 


c. FINITE ELEMENT METHOD LIMITATIONS 
The finite element method was successful for back-propagation of fields measured 


at distances less than 0.25 wavelengths from the surface of the sphere. When the field was 
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measured at distances greater than 0.25 wavelengths, errors due to back-propagation 
increased significantly. The same cases for which the measured field was generated were 
used for back-propagation. These results are summarized in Table 2. The program 
FEM2.M was used to determine the surface current due to back-propagation. It can be 


found in the Appendix. 













BACK-PROPAGATION RESULTS: FINITE ELEMENT METHOD 


Case Sphere Field Modes Jeg 
Radius | Distance RMS Error | RMS Error 
(A) (A) (7%) (%) 


FEM Ila 0.0590 0.9976 
i al 0.0338 30.09 
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Table 2. Centered Sphere Back-Propagation Results: Finite Element Method 






10.26 
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The magnitude and phase of the surface current on the PEC sphere, as determined 
through the finite element method, is shown in Figures 31 through 42. These figures help 
to tell the complete story of the accuracy in back propagating to obtain the current. In 
each figure the “exact” quantity 1s depicted by a solid line and the calculated quantity is 


depicted by dots. Note that in both the magnitude and phase plots the calculated quantity 


closely matches the exact quantity for the cases in which the H field is measured less than 
0.25 wavelengths away from the surface of the sphere. In the cases where the distance is 
one wavelength, the error is immense. Cases FEM2d and FEM3d have been added to the 
original nine test cases in order to demonstrate what happens when the measurement 
distance is increased to 0.26 wavelengths from the surface of the sphere. The following 


section will investigate the cause of this error. 
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Figure 31. Magnitude Comparison: FEM1a 


Jexact (line) & Jback (dots): FEM1a 
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Figure 32. Phase Comparison: FEM1a 
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Jexact (line) & Jback (dots): FEM1b; RMS Error= 30.09% 
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Figure 33. Magnitude Companson: FEM1b 


Jexact (line) & Jback (dots): FEM1b 
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Figure 34. Phase Comparison: FEM1b 
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Jexact (line) & Jback (dots): FEMic; RMS Error= 25.87% 
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Figure 35. Magnitude Comparison: FEM1c 
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Figure 36. Phase Comparison: FEMIc 
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Jexact (line) & Jback (dots): FEM2a; RMS Error= 0.7406% 
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Figure 37. Magnitude Comparison: FEM2a 


Jexact (line) & Jback (dots): FEM2a 
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Figure 38. Phase Comparison: FEM2a 
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Jexact (line) & Jback (dots): FEM2b; RMS Error= 57.18% 
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Figure 39. Magnitude Comparison: FEM2b 


Jexact (line) & Jback (dots): FEM2b 
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Figure 40. Phase Comparison: FEM2b 
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Jexact (line) & Jback (dots): FEM2d; RMS Error= 662.9% 
0.7 


|Jtheta| 
oO Oo 92 2 2 
Nb wo Bh AQ OO 


O 
oo 





2) 


0 90 100 150 200 
Theta (degrees) 


Figure 41. Magnitude Comparison: FEM2d 
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Figure 42. Phase Comparison: FEM2d 
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Jexact (line) & Jback (dots): FEM3a; RMS Error= 0.7333% 
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Figure 43. Magnitude Comparison: FEM3a 
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Figure 44. Phase Comparison: FEM3a 
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Jexact (line) & Jback (dots): FEM3d; RMS Error= 5.547% 
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Figure 45. Magnitude Comparison: FEM3d 
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Figure 46. Phase Comparison: FEM3d 


49 


When back-propagating measurements are made at distances greater than 0.25 
wavelengths from the surface of the sphere using the FEM approach, the RMS error of the 
predicted surface current tends to increase dramatically relative to the situation of smaller 
distances. This appears to result from resonant modes injecting themselves into the 
solution in varying degrees. Figure 47 shows the single mode RMS error for a centered 
PEC sphere used in case FEM1a. Note that the RMS error associated with each mode is 
relatively low; all but two are below one percent. The results shown Table 2 indicate that 
back-propagation gave a solution with less than one percent RMS error. Figure 48 shows 
the single mode RMS error for case FEM1b. In this case the measured field was one 
wavelength from the surface of the sphere. Note that the lower order modes generally 
have RMS error in the neighborhood of five percent or less, except for mode n = 5 which 
is nearly fifty percent RMS error. This error “spike” indicates a resonant mode that 
significantly alters the solution. As can be seen in Figures 48 through 53, the resonant 
modes seem to appear in only the cases where the measurement distance is greater than 


0.25 wavelengths from the surface of the sphere. 
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Single Mode RMS Error for Centered PEC Sphere: FEM2a 
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Figure 49. Single Mode RMS Error: Case FEM2a 


Single Mode RMS Error for Centered PEC Sphere: FEM2b 
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Figure 50. Single Made RMS Error: Case FEM2b 
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Single Mode RMS Error for Centered PEC Sphere: FEM2d 
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Figure 51. Single Mode RMS Error: Case FEM2d 
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Single Mode RMS Error for Centered PEC Sphere: FEM3d 
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Figure 53. Single Mode RMS Error: Case FEM3d 


Another way to view this phenomenon is by looking at the product solution parts 


of H,(r,@) = g(r)-h(@) for a given mode, n. The radial variation, g(r), to be shown in 


the upper subplots, indicates the contribution at the PEC sphere surface due to a constant 
contribution of g(b) = 1 onthe outer boundary. Figures 54 and 55 show the radial 
variation for modes n = 4 and n=5S, respectively, for the case FEM1b. Note that mode 4 
provides relatively small contribution at the PEC sphere surface. Mode 5, however, 
indicates an extremely large contribution at the PEC sphere surface. Two additional cases 
have been added here to further illustrate this observation. Cases FEM2d and FEM3d 
shown in Figures 56 and 57, respectively, are both cases in which the measured field is 
0.26 wavelengths from the surface of the sphere. Note the large contribution at the 


surface of the sphere due to the resonant modes captured in each of these figures. 
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Product Solution Parts of Hphi(r,theta): FEM1b; n=4 
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Figure 54. Product Solution Parts: Case FEM1b, n= 4 


Product Solution Parts of Hphi(r,theta): FEM1b; n=5 
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Figure 55. Product Solution Parts: Case FEMI1b, n= 5 
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Product Solution Parts of Hphi(r,theta): FEM2d; n=8 
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Figure 56. Product Solution Parts: Case FEM2b, n= 8 
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Figure 57. Product Solution Parts: Case FEM3b, n= 16 
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Additional insight is provided by viewing the resonant modes that appear when the 
field-is measured at varying distances from a fixed size sphere. Figure 58 considers modes 


n= 1 to 100. Each line in Figure 58 displays the contribution of g(a) from a single mode 
at the sphere surface as a result of a constant contribution of g(b) = 1 at the boundary 
distance. The “best radius,” 5 = 5.641 A, in this case, provides the minimum g(a) from 
all modes applied as g(b) = 1 . From the perspective of minimal resonance effects, this is 
the best location to measure the field. 


Hphi(a) due to Pn‘“1 at b; a = 5wi; n = 1:100 
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Figure 58. Resonant Effects at Various Outer Boundary Locations 


D. TRANSFER MATRIX FOR CENTERED SPHERE 
As can be seen in Table 2 the results due to back-propagation are poor for 
measurements made greater than approximately 0.25 wavelengths from the surface of the 


sphere. A transfer matrix using spherical harmonics for a centered sphere was developed 


=| 


to accurately back propagate the field on the outer boundary to that on the sphere, 
without the bothersome injection of resonant modes as appear with the FEM solution. 
This method will give some indication as to the spatial resolution required when measuring 
the field so that the source current on the sphere can be accurately determined. 


Figure 58 shows a meridian plane for a centered sphere of radius a enclosed by an 


outer boundary arc of radius 6. For any radius r , where a<r<J, the total H field can 


be determined by: [14] 


Ns Oa: H® (kr 
HG, Oe >» 4, ty an P'(cos@) (37) 
n=] kr kr 
where: 
J; (ka) 
"= FO ay “= 


forces E, = 0 for each mode on the surface of the sphere. 
A specified f(@) = H,(6,@) will provide a set {a,} . The a, coefficients are 


determined by moment matching, using orthogonality of Legendre functions. This 


technique 1s detailed in [14 ]. The transfer matrix H,is developed such that 


H,(,0,) = H,(,P): H,(6,8,) (39) 
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Centered Sphere radius = 2; Arc radius = 5 
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Figure 59. Meridian Plane of Centered Sphere 


where H,(6,0,) is a column array of the @, sampled field at radius 5, as determined in 
Chapter II. The array product results in the column array of H, at radius a which 


corresponds to the surface of the sphere. From this result the surface current is then 
determined using Equation (9b). The test cases using the transfer matrix were generated 
using the program HSPHERE1.M which can be found in the Appendix. Table 3 tabulates 
the results and compares them to those Birained using the finite element metmnod It 
should be noted that the transfer matrix works only for a centered sphere and can not be 


applied to any arbitrary axisymmetric body. 
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SURFACE CURRENT RESULTS USING TRANSFER MATRIX 


Sphere Field Jo 
Radius ]| Distance Segments | RMS Error | RMS Error 
(%) (%) 


la | 1 7. 12 | 4 |  +«»94 ~+&2¥J 0.0590 | 0.8085 0590 0.8085 
0.0338 0.6271 





Table 3. Surface Current Results Using Transfer Matrix 





The magnitude and phase of the surface current on the PEC sphere, as determined 
by the transfer array, is shown in Figures 60 through 77. Once again these figures help to 
tell the complete story of the surface current. In each figure the exact quantity is depicted 
by a solid line and the calculated quantity 1s depicted by dots. Note that in both the 


magnitude and phase plots the calculated quantity closely matches the exact quantity. 
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Jexact (line) & Jback (dots): Caseia; RMS Error= 0.8085% 
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Figure 60. Magnitude Comparison: Casela 


Jexact (line) & Jback (dots): Caseta 
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Figure 61. Phase Comparison: Casela 


|Jtheta| 


Jexact (line) & Jback (dots): Case1b; RMS Error= 0.6271% 
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Figure 62. Magnitude Comparison: Caselb 
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Figure 63. Phase Comparison: Caselb 
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Jexact (line) & Jback (dots): Case1c; RMS Error= 0.7349% 
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Figure 64. Magnitude Comparison: Caselc 
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Figure 65. Phase Comparison: Caselc 


Jexact (line) & Jback (dots): Case2a; RMS Error= 0.7945% 
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Figure 66. Magnitude Comparison: Case2a 
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Figure 67. Phase Comparison: Case2a 
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Jexact (line) & Jback (dots): Case2b; RMS Error= 0.9453% 
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Figure 68. Magnitude Comparison: Case2b 
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Figure 69. Phase Comparison: Case2b 
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Jexact (line) & Jback (dots): Case2c; RMS Error= 0.5991% 
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Figure 70. Magnitude Comparison: Case2c 
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Figure 71. Phase Comparison: Case2c 
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Jexact (line) & Jback (dots): Case3a; RMS Error= 0.7885% 
Omi 


0.1 


© 
oO 
Co 


|Jtheta| 
oO 
oO 
Oo 


0.04 


0 50 400 150 +200 
Theta (degrees) 


Figure 72. Magnitude Comparison: Case3a 
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Figure 73. Phase Comparison: Case3a 
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Jexact (line) & Jback (dots): Case3b; RMS Error= 0.8645% 
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Figure 74. Magnitude Comparison: Case3b 


Jexact (line) & Jback (dots): Case3b 
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Figure 75. Phase Comparison: Case3b 
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Jexact (line) & Jback (dots): Case3c; RMS Error= 0.7846% 
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Figure 76. Magnitude Comparison: Case3c 
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Figure 77. Phase Comparison: Case3c 
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IV. ANALYSIS OF RESULTS 


A. THE INTEGRATION ALGORITHM 

Since the “measured” fields to be back-propagated were to be determined through 
integration, the first part of this thesis was dedicated to generating accurate fields through 
numerical integration. As is shown in Table 1, properly performed numerical integration 
provides highly accurate measured fields. Figures 9 through 12 provide several examples 
of how quickly the integration converges. In each case the percent RMS error is 
determined essentially by the number of points on the sphere surface over which the 
integration takes place. For the cases included in this thesis the number of points was 
chosen such that the error due to back-propagation, as determined by the transfer matrix, 
would be kept below one percent. 

Figures 13 through 30 show the magnitude and phase of the measured field 
determined through numerical integration. For each figure the measured field is over-laid 
on the exact field. Visual inspection indicates the magnitudes are nearly the same over the 
entire range of theta. RMS error computations indicate a difference of only a few 
hundredths of one percent. Visual inspection of the phase comparison gives similar results 
- both the computed and the exact are nearly the same. 

As the number of surface points increase for a constant sphere size, the accuracy 
of the measured H field increases. This result takes the form of decreasing RMS error as 
shown in the integration convergence of Figures 9 through 12. As the measurement 
distance increases for a constant sphere size, the number of points required for accurate 


integration results decreases. This can be seen by comparing Cases 2a and 2b shown in 
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Figures 11 and 12. For the larger measurement distance of Case 2b the RMS error is 


significantly smaller using the same number of points as Case2a. 


B. FINITE ELEMENT METHOD PROBLEMS 

The finite element method was used successfully to determine the source current 
on the PEC sphere when the field was measured at distances less than 0.25 wavelengths 
from the sphere surface. Table 2 summarizes the test cases used in this thesis. 

The large errors that arise when the field to be back-propagated was measured at 
distances greater than 0.25 wavelengths from the sphere surface appear to be caused by 
resonant modes. These modes are shown for several cases in Figures 48, 50, 52 and 53. 
When back-propagated, these modes provide overwhelmingly disproportionate 
contributions at the sphere surface. These contributions are shown for a few 
representative cases in Figures 55, 56, and 57. It appears these resonant modes come 
about as a result of the boundary conditions placed on the outer boundary, the location at 
which the field is measured. As can be seen in Figures 48, 50, 51 and 53, these resonant 


modes appear as sharp spikes in the plot. 


c. TRANSFER MATRIX 

The transfer matrix that was used to determine the source current on the body 
gave results with less than one percent RMS error. The low error is a result of the number 
of surface points used to determine the measured field, and the number of points at which 
the field was measured at the outer boundary. The relationship between number of 
surface points and numerical integration accuracy was discussed in Section A. As can be 


seen in Equation (18), the measured field varies inversely with distance from segment to 
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field point. This inverse distance variation 1s the main point of concern when R is small at 
the closest point of approach. There is a large relative contribution to the integration 
from this surface region and the variation in the field is rapid. The number of segments on 
the sphere surface in this near-region must be increased to ensure the accuracy of 


numerical integration. 
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Vs CONCLUSIONS 


The objective of this research was to investigate the viability of back-propagating 
electromagnetic field measurements to an axisymmetric body in order to determine the 
source distribution for radiated power from the body. The method for achieving this 
objective was straightforward. First, the “measured fields” were simulated through 
numerical integration. Second, these fields were input into a finite element algorithm to 
determine the source current on a PEC axisymmetric body. The original thrust of the 
thesis was to investigate the ranges and resolution at which the fields could be measured in 
order to provide a usable level of accuracy in determining the source distribution on the 
body. The test cases were to include various “arbitrary” axisymmetric shapes to include 
cones, offset cones, offset spheres, and various cone-cylinder-sphere combinations. The 
centered sphere was to be used for very basic test purposes only. 

Along the way to the objective it was found that the finite element method 
provided inaccurate results in the centered sphere tests for all but very limited cases. This 
discovery prompted a redirection of efforts to investigate the sources of error. 
Consequently, progress towards the original goals was sidetracked; but, that is often the 
nature of true research into the unknown. 

The integration program developed to determine the measured fields was found to 
converge to an accurate solution. Although the shape tested was a sphere, the theory 
developed in Chapter II using stacked circular rings can be easily extended to arbitrary 
axisymmetric shapes. 

The finite element method was selected as a means to determine the source current 


on the sphere through back-propagation. The errors associated with this method for the 
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centered sphere were unexpected. Consequently, solving these associated resonance 
problems was not intended to be the focus of this thesis. Further testing of the FEM 
solution will be done as part of extensions to this effort in follow-on thesis efforts and 


techniques for mitigating resonance effects will be explored. 
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APPENDIX - MATLAB CODES 


% thesis2.M by Dan Wawrzyniak 5/10/97 

% updated 5/14/97 

% Combines all programs into one 

% Includes checking integration error first 

% Allows user to save data in a .mat file for future use 
clear all; 

case = input(‘What case are you running? ','s’); 

fO0 =1input(Enter frequency in MHz: '); 

lambda = 3e2/f0; 

disp({'Wavelength = 'num2str(lambda),' meters']); 
al =input(‘Enter metal sphere radius in terms of wavelength: ');a = al*lambda; 
rl =input(Enter outer boundary in terms of wavelength: ');r =r1*lambda; 
eta =120*pi; % free-space Z_0 

Idl =1; % dipole moment 

B =2*pi/lambda; % wavenumber 

z01 = input(Enter z-position of elemental dipole in terms of wavelength: '); 
zO = z01*lambda; 

Nth1 = input(Enter number of theta points: ‘); Nth = Nth1-1; 

ppwl = input(‘Enter number points per wavelength: '); 

dth=pi/Nth; th=(0:dth:pi)'; thdeg=180*th/pi1; 

sth=sin(th); cth=cos(th); 

% computing rd, sin(thd), and cos(thd) for sphere surface (r=a) 
rd=sqrt(a*at+z0*z0- 2*a*z0*cth); sthd=a*sth./rd; cthd=(a*cth-z0)./rd; 

% computing field components in dipole centered Se 
Fth=Id|*exp(-}*B*rd)./(4*pi*rd); 

Hpi=Fth.*sthd.*(qj*B + 1./rd); 

Etd=eta*Fth. *sthd.*(j*B + 1./rd - j./(B*rd.*rd)); 
Erd=2*eta*Fth.*cthd.*(1./rd - j./(B*rd.*rd)); 

% translating incident field spherical components 

cdth=cth.*cthd + sth.*sthd; sdth=sth.*cthd - cth.*sthd; 

Eti=Etd.*cdth - Erd.*sdth; % note that Hpi=Hpd 

% Computing spherical harmonic coefficients for scattered field 

% Ets expansion which minimizes the total sphere surface 

% tangential field = Etit+Ets 

Nmax = fix(2*B*r); 

disp({"Nmax = ';num2str(Nmax)]}); 

N = input(Enter number of spherical harmonics: '); 

Sn=-dth*Eti.*sth; Pnl=legpol2(cth,N); In=(Pn1.')*Sn; 

n=(1:N)'; cn=(2* n+1)./(2*n.*(n+1)); 

Ba=B*a; [Hn,DHn]=shan3(Ba,N); 

Hn = Hn.'; DHn = DHn.'; 

Dn=DHn/Ba; an=(cn.*In)./(j*eta*Dn); 

Ets=j*eta*Pn1 *(an.*DHn)/Ba; Ett=EtitEts; % ideally Ett --> 0 
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Hps=Pn1 *(an.*Hn)/Ba; Hpt=Hpit+Hps; % Jtheta=H_pt 
Jtheta = -Hpt; 
Br=B*r; [Hn, DHn]=shan3(Br,N); 
Hn = Hn.'; DHn = DHn.'; 
Ets=j*eta*Pn1*(an.*DHn)/Br; Hps=Pn1*(an.*Hn)/Br; 
numring = Nth; % number of rings 
inc = input('Enter number of face segments: '); 
ppoints = input(‘Enter phipoints - wavelength factor: '); 
% Field Point Inputs 
% Input the Field Point in terms of theta. 
radius =r; % radius of field point 
numFP = length(Hps); 
thetaFP = linspace(0,pi,numFP); 
thetaFPd = 180*thetaFP/p1; 
rho = radius. *sin(thetaFP); 
Z = radius. *cos(thetaFP); 
% Jt is the tangential surface current at the nodes. 
% The nodes are defined by rhop & zp. 
Jt = Jtheta; 
% Determine the thetas for the sphere 
dtheta = pi/numring; % gives the delta theta 
theta = O:dtheta:pi; % gives numring +1 values of theta. 
thetad = theta.*180/pi; % theta in degrees 
“% Determine (rhoprime,zprime) pairs on the sphere 
rhop = a*sin(theta); 
Zp =a*cos(theta); 
% Find the length "s" of each face. 
s = sqrt(rhop(2)*2 + (zp(1)-zp(2))*2); 
% Find the angle associated with each face 
alpha = acos(rhop(2)/s) + theta(1:numring); % radians 
alphad = alpha* 180/pi1; % degrees 
talpha = tan(alpha); 
alphal = ones(inc+1,1)*alpha; | 
alpha = reshape(alphal ,numring*(inc+1), 1); 
dzf = (zp(1:numring)-zp(2:numring+1 ))/inc; 
phipoints = round(p1.*a*ppoints/lambda)+1; 
dphi =pi./phipoints; 
phi = 0:dphi:p1; 
cp = cos(phi); 
for n= 1:numring; % loops through rings 
zf(:,n) = (zp(n):-dzf(n):zp(n+1)); 
rhof(:,n) = (linspace(rhop(n),rhop(n+1),inc+1))’; 
JT¢:,n) = Jt(n)*(zfC.:inc+1,n)-zp(n+1))/(zp(n)-zp(n+1))+-... 
Jt(nt+1)*(zf( :inc+1 ,n)-zp(n))/(zp(n+1)-zp(n)); 
end 
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% Reshape the matrices to column vectors 
zf = reshape(zf,numring*(inc+1), 1); 
rhof = reshape(rhof,numring* (inc+1),1); 
JT =reshape(JT,numring*(inc+1), 1); 
for FP = 1:numFP; % loops through field points 
Rl = 2*rho(FP)*rhof*cp; 
R2 = (rho(FP)*2+rhof.*2+(z(FP )-zf).%2)*ones(1,length(cp)); 
R = sqrt(R2-R1); 
[Fcs,F1s] = fcesfls2(cp,dphi,R); 
Ha = rhof.*sin(alpha); 
Hb = (z(FP) - zf).*cos(alpha); 
Hc = (Ha-Hb).*Fcs; 
Hd = rho(FP).* sin(alpha).*F 1s; 
H =(Hc-Hd).*JT; 
H = reshape(H,inc+1,numring); 
Fil = a*((1/(4*pi))*(trapz(H))); 


H1 = dzf.*H1; 
HH(1:numning,FP) = H1.'; 
end 
Hphi = sum(HH); 


%% RMS Error Calculations 

E = 100*sqrt((Hps - Hphi.')*(Hps - Hphi.'))./sqrt(Hps'* Hps); 

figure 

plot(thdeg,abs(Hphi),'ko',thdeg,abs(Hps),'k') 

title(["Hps (exact: -, integrated: 0): ',eval(‘case’),'; RMS Error= ' ;num2str(E),'%']) 

xlabel('Theta (degrees)'); 

ylabel(‘|Hscattered|'); 

figure 

plot(thdeg,real(Hphi),'k*',thdeg,imag(Hphi),'k+'’,thdeg,real(Hps).... 
'k',thdeg,imag(Hps),'ko’); 

title({"Hexact(real: -, imag: o) & Hint(real: *, imag: +): ',eval(‘case')]}) 

xlabel(‘Theta (degrees)') 

ylabel(‘Hps (Amps/meter)') 

Ytext(.55,.95,['Frequency= ';num2str(f0),'Mhz’'],'color’,... 

% 'g','units',"normalized’); 

“%text(.55,.90,['Wavelength= 'num2str(lambda),'meters'],'color',... 

% 'g'units',{normalized’), 

“%text(.55,.85,['Dipole Offset= ';num2str(z01),'*lambda'],'color’,... 

% ‘g',units’,{normalized’), 

Ytext(.55,.80,['Radius of Sphere= ',;num2str(al),'*lambda'],'color',... 

%  'g',units',/normalized'); 

%text(.55,.75,['Radius of Field Point= ';num2str(r1),'*lambda'],'color'.... 

% 'g'units,,'normalized'); 

“%text(.55,.70,["Number Rings= ',;num2str(numring)],'color'".... 

%  'g',‘units','normalized’); 


81 


%text(.55,.65,["Face Segments= ',num2str(inc)],'color’,... 
% 'g'‘units','normalized’); | 
%text(.55,.60,["Phi-Wavelength Factor= ';num2str(ppoints)],'color',... 
% 'g' ‘units','normalized'); 
%text(.55,.55,["Field Points= ',;num2str(numFP)],'color’,... 
% 'g'‘units',‘normalized’); 
Ytext(.55,.50,['Field Points per wl= ',num2str(ppwl)],'color',... 
% 'g'‘units',‘normalized’'); 
“text(.55,.45,["Number modes= ';num2str(N)],'color",... 
% '‘g'‘units',‘normalized'); 
%text(.55,.40,['RMS Error= ';num2str(E),' %'],'color',... 
%  'r',‘units','normalized'); 
M1100 %%%0% 
% Computes Hphi incident at the field point radius due to the dipole. 
% This result is added to the Hphi scattered computed using Dan's 
% integration program 
dth=pi/Nth; th=(0:dth:pi)'; thdeg=180*th/p1; 
sth=sin(th); cth=cos(th); 
% computing rd, sin(thd), and cos(thd) for sphere surface (r=a) 
rd=sqrt(r*r+z0*z0-2*r*z0*cth); sthd=r*sth./rd; cthd=(r*cth-z0)./rd; 
% computing field components in dipole centered coordinates 
Fth=Idl* exp(-j* B*rd)./(4*pi*rd); 
Hpi=Fth.*sthd.*(4*B + 1./rd); 
Hphi = Hphi.'; % make Hphi a column vector 
MAHWAH’ AWAWDVWAV1%%%%0%%%%%%%%%%%%%%%%%%% 
% This program loads the transfer function Hs from femdat.mat 
% and created by FEM2.M and dots it with Hphi stored in 
% source.mat from Dan's integration program SOURCE.M. 
% 
% The result is the current on the surface of the sphere 
% arrived at through the back propagation of the scattered 
% H field. 
% 
% The result (Jback) is compared to the surface current used in 
% Dan's program (Jtheta, stored in jtheta.mat). 
% 
% The surface current is generated by DIPSPHR2.M 
% 
load femdat.mat 
“ Add Hphi scattered from sphere and Hpi incident from dipole 
Hphitot = -(Hphi + Hpi); 
“% Determine Surface current due to back-propagation of Hphi 
% Hhpi is UNFILTERED 
Jback = Hs*Hphitot; 
“% Determine Surface current due to back-propagation of Hphi 
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% RMS Error Calculations 
E = 100*sqrt((Jtheta - Joack)'*(Jtheta - Jback))./sqrt(Jtheta'*Jtheta); 
figure 
plot(thdeg,abs(Jtheta),'k',thdeg, abs(Jback),'ko') 
title(['Jexact (-) & Jback (0): ',eval(‘case'),'; RMS Error‘... 
num2str(E),'%']) 
xlabel(‘Theta (degrees)') 
ylabel('|Jtheta]") 
%legend(‘Jtheta - dipsphr2.m','Jback: no filtering - fem2.m') 
figure 
plot(thdeg,real(Jtheta),’k',thdeg,imag(Jtheta),'ko',thdeg,real(Jback).... 
'k*' thdeg,imag(Jback),'k+') 
title([‘Jexact(real: -, imag: o) & Jback(real: *, imag: +): 'eval(‘case')]); 
xlabel(‘Theta (degrees)') 
ylabel(‘Jtheta (Amps/m*m)') 
yn = input(‘Save data from this run? (Y/N): ','s'); 
if ~isempty(yn) 
if yn —_— 'y' yn —_— iy 
clear eta Idl B a z0 dth th sth cth rd sthd cthd Fth Etd Erd cdth sdth ... 
Eti r Sn Pnl Inn cn Ba Hn DHn Ets Br Dn Ett Fls FP Fcs E 
clear H H1 HH Ha Hb Hc Hd JT Jt R R1 R2 alpha alphal alphad an cp dphi... 
dtheta dzf rb fho rhof rhop rs s talpha yn z zb zf zp zs 
clear Hpt Nmax Nth f phi phipoints... 
radius rho theta thetaFP thetaFPd thetad 
save (eval('case')) 
%save f0 lambda al z01 rl Nth N inc ppoints Hs Hps numning ppoints numFP... 
%  ppwl name thdeg Hps Jback Hphi Hphitot 
end 
end 
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function [Fc_s,Fl_s] = fcsfls2(cp,dphi,R) 
% For use with HPHIERR2.M 

% 

% This function does the integration to find Fc(s) & F1(s) 
% updated 4/19/79 at home Dan Wawrzyniak 
lambda = 1; 

beta = 2* pi/lambda; 

phase = exp(-i*beta.*R); 

same = (1+i*beta.*R)./R.43; 

F1 = same. * phase; 

Fl_s = (2*dphi*(trapz(F1.'))).'; 

Fc = ones(size(F1_s))*cp.*same.* phase; 

Fc_s = (2*dphi*(trapz(Fc.'))).'; 
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% DipSpr2.M computes scattered fields due to an elemental z-axis 
% dipole positioned above the north pole of a metallic sphere. 

% Mod-2 by M.A. Morgan 3/14/97 

% shan3.m & legpol2.m loaded 5/2/97 


clear all; 
eta=120*pi; % free-space Z_0 
Idj=1; % dipole moment 


f0=300; % 300 MHz 
B=2*pi*f0/300; % wavenumber 
z0=input(‘Enter z-position of elemental dipole in meters: '); 
a=input(‘Enter metal sphere radius in meters: '); 
Nth=input(Enter number of theta segments: '); 
dth=pi/Nth; th=(0:dth:pi)'; thdeg=180*th/pi; 
sth=sin(th); cth=cos(th); 
% computing rd, sin(thd), and cos(thd) for sphere surface (=a) 
rd=sqrt(a*a+z0*z0-2*a*z0*cth); sthd=a* sth./rd; cthd=(a*cth-z0)./rd; 
% computing field components in dipole centered coordinates 
Fth=Idl*exp(-j*B*rd)./(4* pi*rd); 
Hpi=Fth.*sthd.*(j*B + 1./rd); 
Etd=eta*Fth.*sthd.*(j*B + 1./rd - j./(B*rd.*rd)); 
Erd=2*eta*Fth. *cthd.*(1./rd - j./(B*rd.*rd)); 
% translating incident field spherical components 
cdth=cth.*cthd + sth.*sthd; sdth=sth.*cthd - cth.*sthd; 
Eti=Etd.*cdth - Erd.*sdth; % note that Hpi=Hpd 
% Computing spherical harmonic coefficients for scattered field 
% Ets expansion which minimizes the total sphere surface 
% tangential field = Etit+Ets 

=input(‘Enter number of spherical harmonics: '); 
Sn=-dth*Eti.*sth; Pnl=legpol2(cth,N); In=(Pn1.')*Sn; 
n=(1:N)'; cn=(2*n+1)./(2*n.*(n+1)); 
Ba=B*a; [Hn,DHn]=shan3(Ba,N); 
Hn = Hn.'; DHn = DHn..'; 
Dn=DHn/Ba; an=(cn.*In)./G*eta*Dn); - 
Ets=j*eta*Pn1 *(an.*DHn)/Ba; Ett=Etit+Ets; % ideally Ett --> 0 
Hps=Pn1*(an.*Hn)/Ba; Hpt=Hpi+Hps; % Jtheta=H_pt 
Jtheta = -Hpt; 
save dipsphr2 Jtheta thdeg 
% Plotting expansions on sphere surface 
clf reset; plot(thdeg,abs(Eti),thdeg,abs(Ett),’.'); 
xlabel(‘theta (deg)'); 
ylabel(‘|E_thetal: Incident (solid); Total (dots)’); 
title({‘Dipole at zO=';num2str(z0),'; Sphere Radius=',num2str(a).... 
'; No. Modes="'int2str(N)]); 
figure(1); 
hcpy=input(Enter 1 for Hard Copy: '); 
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if hcpy == 1, print; end; 

clf reset; plot(thdeg,abs(Hpi),thdeg,abs(Hps),’.',thdeg,abs(Hpt),'--'); 
xlabel(‘theta (deg)'); 

ylabel(‘|H_ phil: Incident (solid); Scattered (dots); Total (dashed)’); 
title({‘Dipole at zO=',num2str(z0),'; Sphere Radius=',num2str(a).... 
' No. Modes="',int2str(N)]); 

figure(1); 

hcpy=input(‘Enter 1 for Hard Copy: '); 

if hcpy = 1, print; end; 

% Now compute scattered fields at specified radius 
r=input(‘Enter radius in meters to field point: '); 

Br=B*r; [Hn,DHn]=shan3(Br,N); 

Hn = Hn.'; DHn = DHn.'; 

Ets=j*eta*Pn1 *(an.*DHn)/Br; Hps=Pn1 *(an.*Hn)/Br; 

% Plotting scattered field expansions 

clf reset; plot(thdeg,abs(Ets),thdeg,eta*abs(Hps),'.'); 

xlabel(‘theta (deg)’); 

ylabel(‘|E_ theta] (solid); eta*/H_phij (dots)’); 

title({"Scattered Fields at r=',num2str(r),' m; Dipole at z0=",... 
num2str(z0),'; Sphere Radius=",num2str(a),'; No. Modes="',int2str(N)]); 
figure(1); 

hcpy=input('Enter 1 for Hard Copy: '); 

if hcpy = 1, print; end; 

end; 
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% FEM2.M provides finite-element solution for user-specified 
% metallic body of revolution using VARELA algorithm. 

% By M.A. Morgan 3/20/97- 4/17/97 

clear all 

% Generating surface and boundary points then mesh topology 
sgen=input(‘Enter name of surface generation program: ','s'); 
eval(['[rs zs rb zb]=',sgen,';']); 

{rho zee nd elnd]=mesh2(rs,zs,rb,zb); % Computing mesh 

% Computing Mesh Parameters 


N=length(rho); % Number of Nodes in Mesh 
L=length(elnd(:,1)); | % Number of Elements in Mesh 
Ns=length(nd); % No. of Surface or Boundary Nodes 


Nm=N-nd(1)-nd(Ns)-Ns+2; % Number of Nodes for Unknown Hphi 
disp(’’) 
disp(‘Finite Element Mesh Parameters: ') 
disp({’ Number of Nodes: ‘int2str(N)]) 
disp({' Number of Elements: 'int2str(L)]) 
disp({' Number of Unknowns: ',int2str(Nm)]) 
disp(’ ’) 
disp(‘Press a Key to Display Mesh or <Ctl> C to Abort’); 
pause 
% Plotting Mesh Nodes and Elements 
“oclf reset 
figure('PaperPosition',[1.5 0.5 5 3.75]) 
if Ns <=19, plot(rho,zee,'go’); 
else, plot(rho,zee,'k.'); end 
hold on 
for 1=1:L 
nds=elnd(1,:); 
rl=[rho(nds); rho(nds(1))]; zl=[zee(nds); zee(nds(1))]; 
plot(rl,zl,'k’); 
end 
v=axis; V(1)=0; v(2)=v(4)-v(3); axis(v); axis square; 
title(‘2-D Axisymmetric Shape with mesh’) 
xlabel('meters') 
ylabel(‘meters’) 
text(2.7*v(2)/5,5*v(4)/6,'Mesh Parameters: ') 
text(3 *v(2)/5,3 *v(4)/4,['Nodes: ',int2str(N)]) 
text(3 *v(2)/5,2*v(4)/3 ,['Elements: 'int2str(L)]) 
text(3*v(2)/5,7*v(4)/12,['Unknowns: ',int2str(Nm)]) 
yn=input('Print Hard Copy ? (Y/N): ','s'); 
if ~isempty(yn), if yn == 'y' | yn == 'Y', print; end; end 
% Loading System Matrices By Indexing Through Elements 
disp(‘Sparse Matrix Loading Using Element Contributions Can Now Begin’) 
f=input(‘Enter Frequency in MHz to Start Loading; <Ctl> C to Stop: '); 
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k2=(pi*f/150)*2; T=zeros(3,3); 
=sparse(Nm,Nm); B=sparse(Nm,Ns-2); 
12=0; nf=cumsum(nd); ns=nf-nd+1; % node finish/start per row 
for 1=2:Ns; 
11=12+1; 12=11+nd(i)+nd(i-1)-3; % start/finish elements 
for 1=11:12 
nds=elnd(I,:); rp=rho(nds); zp=zee(nds); 
T=varela(k2,rp,zp); 
disp({"Loading Element ',int2str(1),' of ',int2str(L)]) 
% Loading sparse matrix contributions 
for m=1:3; 
ma=0; 
if nds(m) ~= nf(i-1) & nds(m) ~=nf(i), % m not on boundary 
if nds(m) > nf(i-1) & i< Ns, 
ma=nds(m)-ns(2)-1+3; end 
if nds(m) < ns(i) &1> 2, 
ma=nds(m)-ns(2)-i+4; end 
if ma > 0, % mis a solution node 
for n=1:3; 
na=0; 
if nds(n) ~= nf(i-1) & nds(n) ~=nf(), 
if nds(n) > nf(i-1) & i< Ns, 
na=nds(n)-ns(2)-1+3; end 
if nds(n) <ns(i) &1> 2, 
na=nds(n)-ns(2)-1+4; end 
if na > 0, 
A(ma,na)= A(ma,na)+T(m,n); end 
end % n not on boundary end 
if nds(n) = nf(i-1) & 1> 2, 
B(ma,i-2)=B(ma,i-2)-T(m,n); end 
if nds(n) == nf(i) &1<Ns, 
B(ma,i-1)=B(ma,i-1)-T(m,n); end 
end %n-loop end 
end %ma>Oend 
end % m not on boundary end 
end % m-loop end 
end % |-loop end 
end % i-loop end 
disp(‘Matrices Loaded, Now Solving System ... Be Patient’) 
“% Solve sparse matrix system to construct transfer function 
Yo atray relating internal Hphi to Ns-2 nodal boundary values 
H=A\B;%clear ABCD % freeing memory 
“ Extract rows to form transfer array relating Ns nodal values 
% of Hphi (including two z-axis Hphi=0 BC's) to Ns nodal values 
% of Hphi on the PEC surface. Note that the PEC surface output 
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% column vector has Hphi(1)=Hphi(Ns)=0 since these correspond 
% to nodes on the z-axis. 

Hs=zeros(Ns,Ns); 

row=1; Hs(2,2:Ns-1) = H(row,:); 

for n=3:Ns-1; 

row=row+tnd(n-1)-1; Hs(n,2:Ns-1) = H(row,:); 

end 
n=1:Ns; clf; plot(n,diag(Hs),'r') 
v=axis; axis([1 Ns v(3) v(4)]) 
xlabel(‘row(n)'); ylabel(‘Real diag[Hs]') 
title(["Diag[Hs] for ',sgen,'; Ns=',int2str(Ns),'; N=',int2str(N).... 

'' Nm="',int2str(Nm),'; L="',int2str(L)]); figure(1) 

yn=input(‘Print Hard Copy ? (Y/N): 5,'s'); 
if ~isempty(yn), if yn = 'y' | yn = 'Y'’, print; end; end 
% Saving Needed Parameters and Hs Transfer Function Array 
save femdat f rs zs rb zb Hs 
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% Program FEMCHK4.M compares exact and FEM2.M computed 
% nodal values for spherical harmonic single mode solutions 
% of an offset metal sphere. By M.A. Morgan 4/21/97 
% Mod-2a (4/30/97) uses MatLab Bessel and Hankel Functions in shan3.m 
% Mod-3 (5/2/97) uses Wronskian to simplify r=a "exact" solution 
% and computes and plots errors for applied modes 
% Mod-4 (5/3/97) uses real mode BC's Pn“1(cos(theta)) on r=b 
% which produces a real internal solution 
clear all; 
case = input("What case are you running? ','s’); 
load femdat; % loading frs zs rb zb Hs 
k=pi*f/150; Ns=length(rs); 
d=(zs(1)+zs(Ns))/2 - (zb(1)+zb(Ns))/2; % computing offset 
a=zs(1)-d; b=zb(1); Ra=k*a; Rb=k*b; “% reset sphere centers 
zbd=zb-d; r_b=sqrt(rb.42+zbd.%2); b=zb(1); 
th=linspace(0,pi,Ns)'; c_a=cos(th); th_a=linspace(0,180,Ns)'; 
c_b=zbd./r_b; th_b=180*acos(c_b)/p1; 

=input(‘Enter Upper Spherical Harmonic Mode to Check: '); 
{hna, Dhna}j=shan3(Ra,N); [hnb, Dhnbj=shan3(Rb,N); 
DJna=real(Dhna); Jnb=real(hnb); 
DNna=-imag(Dhna); Nnb=-imag(hnb); 
Dn=(Jnb.*DNna-DJna.*Nnb); An=ones(Ns, 1)*DNna; Bn=ones(Ns, 1 )*DJna; 
(Hn,DHn]=shan3(k*r_b,N); Jn=real(Hn); Nn=-imag(Hn); 
pna=legpol2(c_a,N); pnb=legpol2(c_b,N); pct=zeros(N, 1); 
Ha=b* pna./(a* ones(Ns, 1 )*Dn); % "exact" r=a 
Hb=b* pnb. *(An.*Jn-Bn.*Nn)./(r_b*Dn); % "exact" r=b 
Hc=Hs*Hb; % Computed solution on r=a 
for n=1:N; 
pcet(n)=100* sqrt((Ha(:,n)-He(:,n))'*(Ha(:,n)-Hc(:,n)))/... 

sqrt(Ha(:,n)'*Ha(:,n)); 

end 


[tseg] ,tseg2] = size(Hs); 
figure('PaperPosition',[1.5 0.5 5 3.75]) 
bar(1:N,pct,'k'), v=axis; v(1)=0; v(2)=N+1; v(3)=0; axis(v); 
xlabel(‘Mode Index "n"'); ylabel(‘Percent’); 
title({'Single Mode RMS Error for Centered PEC Sphere: ',eval(‘case’)]) 
yn=input('Print Hard Copy ? (Y/N): ','s'); 
if ~isempty(yn), if yn == 'y' | yn == 'Y', print; end; end 
while 1, 
=input(["Enter n <= ‘int2str(N).... 
‘for Spherical Harmonic (0 to stop): ‘]); 
if isempty(n), break; 
elseif n == 0, break; 
elseif n > N, disp('n exceeds N'); break; end 
figure(‘PaperPosition',[1.5 0.5 5 3.75]) 
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plot(th_b,Hb(:,n),’k') 
xlabel(‘Theta (deg)'); ylabel(“Hphi (A/m)'); 


title((‘Single Mode Hphi at Outer Boundary: ',eval(‘case’).... 


'. n="',int2str(n)]) 
v=axis; axis({O 180 v(3) v(4)]); 
yn=input(‘Print Hard Copy ? (Y/N): ','s'); 
if ~isempty(yn), if yn = 'y' | yn == 'Y', print; end; end 
if length(th_a) < 38, 
figure(‘PaperPosition',[1.5 0.5 5 3.75]) 
plot(th_a,Ha(:,n),'k',th_a,Hc(:,n),'k.') 
title(["Exact (line), FEM (dots): ',eval(‘case’).... 
' n=',int2str(n),'; RMS Error=',num2str(pct(n)),'%'}) 
else, 
figure(‘PaperPosition',[1.5 0.5 5 3.75]) 
plot(th_a,Ha(:,n),'k',th_a,Hc(:,n),'k.'); 
title({"Exact (line), FEM (dots): '‘,eval(‘case’).... 
' n=',1nt2str(n),'; RMS Error=',num2str(pct(n)),'%']) 
end 
xlabel(‘Theta (deg)'); ylabel(‘Hphi (A/m)’'); 
v=axis; axis({O 180 v(3) v(4)]); 
=input(‘Print Hard Copy ? (Y/N): ','s'); 
if ~isempty(yn), if yn = 'y' | yn == 'Y', print; end; end 
end 


91 


% hspherel.m 
% Program Hsphere.M uses spherical harmonics to compute 
% the Hs array for a centered PEC sphere. Saves in same 
% format as FEM2.m and can be used in FEMCHKxx.m programs 
% By M.A. Morgan 5/2/97 
% 5/10/97 modified by Dan Wawrzyniak 
% Nmax = 2*k*b 
clear all 
f=input('Enter frequency in MHz: '); k=pi*f/150; lambda = 3e2/f; 
a=input(‘Enter metal sphere radius "a" in meters: '); 
b=input(‘Enter outer mesh radius "b" in meters: '); 
Ns1=input(‘Enter number of theta points per wavelength: '); 
Ns = fix(pi*b*Ns1/lambda); 
disp({'Total theta points = ',int2str(Ns)]) 
th=linspace(0,pi,Ns); thd=linspace(0,180,Ns); dth12=pi/(12*(Ns-1)); 
Ra=k*a; Rb=k*b; cth=cos(th); sth=sin(th); Nmax=2*fix(Rb); 
disp(' '); disp({'Estimated Nmax=integer[2*k*b] 1s ',int2str(Nmax)]) 
=input(‘Enter Upper Mode Order to Use in Computing Hs: '); 

disp(' ') 
{Hna,DHna]=shan3(Ra,N); [Hnb, DHnb]=shan3(Rb,N); 
Pn1=legpol2(cth,N); n=1:N; Cn=(2*n+1)./(2*n.*(n+1)); 
Qn=(b/a)./(real(DHna).*imag(Hnb)-real(Hnb).*imag(DHna)); 
n=1:N; A=zeros(N,Ns); Hs=zeros(Ns,Ns); 
% Numerical integration of up(th)*Pn1(cth)*sth 
for p=2:Ns-1; 
fn=Pn1(p-1:p+1,:); g=sth(p-1:p+1); 
Inp=dth12*(fn(1,:)*g(1)+6*fn(2, :)*2(2)+n(3,:)*g(3)+... 

fn(1,:)*g(2)+fn(2,:)* g(1)+fn(2,:)*g(3)+in(GG, :)* g(2)); 
AC.,p)=(Cn.*Qn.*Inp).'; 
end 
Hs=Pn1*A; 
n=1:Ns; clf, plot(n,diag(Hs),'r') 
v=axis; axis({1 Ns v(3) v(4)]) 
xlabel(‘row(n)'); ylabel(‘Real Hs Diagonal") 
title({'Diag[Hs] using Hsphere.m for Ns=',int2str(Ns).... 

'; Nmax="',int2str(N), }); 

figure(1) 
yn=input('Print Hard Copy ? (Y/N): ','s’); 
if ~isempty(yn), if yn == 'y' | yn = 'Y', print; end; end 
rs=a* sth’; zs=a*cth'; rb=b*sth'; zb=b*cth’'; 
% Saving Needed Parameters and Hs Transfer Function Array 
save femdat frs zs rb zb Hs 


a 


function [Hn,DHn] = shan3(X,N) 

% Riccati Spherical Hankel Function Hn*(2)(X) = Jn(X) -j* Yn(X) 
% and Derivative. Using MatLab cylindrical functions of order n+1/2. 
% Function returns: Hn(n)=Hn“(2)(x) and DHn(n)=dHn“(2)(x)/dx 
% for n=l to N. Input N is a positive integer. X is an array 

% of real or complex values [x1 x2 ... xM]. Output Hn(xm) and 
% DHn(xm) are M by N complex arrays. By M.A. Morgan. 

% Mod-3 (5/1/97) allows X to be either row or column array. 

[m1 m2]=size(X); M=max(m1,m2); Hn=zeros(M,N); DHn=Hn; 

if m1==1, x=X'; else, x=X; end 

Nu=(1:N)+0.5; xn=x*ones(1,N); 

JO=sqrt(pi*x/2). *besselj(0.5,x); 

Jn=sqrt(pi*xn/2).*besselj(Nu,x); 

N0=saqrt(pi*x/2). *bessely(0.5,x); 
Nn=sqrt(pi*xn/2).*bessely(Nu,x); 

Hn=Jn-j* Nn; HO=J0-j*NO; DHn(:, 1)=HO-Hn(.,1)./x; 

if N > 1, DHn(:,2)=Hn(., 1)-2*Hn(:,2)./x; end 

if N > 2, 

for n=3:N; DHn(:,n)=Hn(.,n-1)-n*Hn(.,n)./x; end 

end 
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function Pn1 = legpol2(X,N) 

% Computing a matrix of associated Legendre polynomials 

% Pn“m(x) with m=1, for n=1, 2, .. N where N is a positive integer 
% and X is an M-element array of real values [x1 x2 ... xM]. 

% Computed real array Pnl(xm) has M rows and N columns. 

% Using upward recurrence formula from page 953 of Balanis - 

% Advanced EM Engineering, Wiley, 1989. Program by M.A. Morgan 
% Mod-2 (5/1/97) allows either row or column input X array 

[m1 m2]=size(X); M=max(m1,m2); Pnl=zeros(M,N); 

if ml=—1, x=X'; else, x=X; end 

Pnl=zeros(M,N); Pnl1(:, 1)=-sqrt(1-x.*x); 

if N > 1, Pn1(.,2)=3*x.*Pnl(:,1); end 

if N > 2, 

for n=2:N-1; 

nl=1/n; Pn1(:,n+1)=(2+n1)*x.*Pnl(:,n)-(1+n1)*Pn1(:,n-1); 

end 

end 
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function [rs,zs,rb,zb] = osphere 

% Computing meridian surface coordinates for offset sphere 
% (rs,zs) with a=radius and d=z-axis offset. Outer 

% mesh boundary is centered b=radius sphere (rb,zb). 

% Returning coordinates in Ns x 4 array. 

% By M.A. Morgan 3/24/97 

disp(‘Offset Sphere Surface and Mesh Boundary Program’); 
a=input(‘Enter sphere radius (meters): '); 

d=input(Enter sphere offset (meters): '); 

b=input(‘Enter mesh radius (meters): '); 

Ns=input(‘Enter number of surface points: '); 
tha=linspace(0,pi,Ns)'; sta=sin(tha); cta=cos(tha); 

rs=a* sta; zs=a*ctatd; 

thb=tha-asin(d*sin(pi-tha)/b); stb=sin(thb); ctb=cos(thb); 
rb=b* stb; zb=b*ctb; 

clf reset; plot(rs,zs,rs,zs,'.g',rb,zb,rb,zb,'.g') 

v=axis; v(1)=0; v(2)=v(4)-v(3); axis(v); axis square; 
xlabel(‘Press a Key to Display Mesh or <Ctrl> C to Abort') 
figure(1) 
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function [rs,zs,rb,zb] = cone 
% Computing meridian surface coordinates for centered cone 
% (rs,zsS) with h=height and b=base radius . Outer 
% mesh boundary is centered a=radius sphere (rb,zb). 
% Returning coordinates in Ns x 4 array. 
% By M.A. Morgan 3/27/97 
disp(‘Conical Surface and Spherical Mesh Boundary Program’); 
h=input(‘Enter cone height (meters): '); 

=input(‘Enter cone base radius (meters): '); 
a=input(‘Enter outer mesh spherical radius (meters): '); 
Ns=input(‘Enter number of surface points: '); 
rs=zeros(Ns, 1); zs=rs; th=rs; 
S=b+sqrt(b*bt+h*h); % total surface length on cone 
Nb=fix(Ns*b/S)+1; 1f Nb>Ns-2, Nb=Ns-2; end; db=b/Nb; Nz=Ns-Nb; 
rs(1:Nz)=linspace(0,b,Nz); zs(1:Nz)=linspace(h/2,-h/2,Nz); 
rs(Nz+1:Ns)=linspace(b-db,0,Nb); zs(Nz+1:Ns)=(-h/2)* ones(Nb, 1); 
% distort boundary node positions to improve mesh near corner 

=atan(h/b)/2; gm=a2+pi/2; c=(h/2)-b*tan(a2); 

di=asin(c*sin(gm)/a); the=gm+dl; dtb=(pi-thc)/Nb; 
th(1:Nz)=linspace(0,thc,Nz); th(Nz+1:Ns)=linspace(thc+dtb, pi, Nb); 
st=sin(th); ct=cos(th); rb=a* st; zb=a* ct; 
clf reset; plot(rs,zs,rs,zs,'r.',rb,zb,rb,zb,'g.') 
v=axis; v(1)=0; v(2)=v(4)-v(3); axis(v); axis square; figure(1) 
xlabel('Press a Key to Display Mesh or <Ctrl> C to Abort') 
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